#!/usr/bin/env python
# -*- coding: iso-8859-1 -*-
# This module is mostly of untested work, that tries out several modules
#already implemented by Sage. Some are possible of being consolidated to curves.py
from sage.all import *

# Works only for curves with h polynomial equals zero.
def curveWithoutSingularPoints(f):
    return f.derivative().gcd(f) == 1


def plotCurve(fstr, precisao = 350, xRange = 2, yRange = 3):
    x = polygen(QQ)
    f = eval(fstr)
    h = 0 * x
    func = lambda px, py: py ** 2 + h(px) * py - f(px)
    c = contour_plot(func, (-xRange, xRange), (-yRange, yRange),
                     plot_points = precisao, cmap = 'winter', fill = False, contours = [0])
    R = line([(0, -yRange), (0, yRange)], linestyle = '--', rgbcolor = (0, 0, 0))
    L = line([(-xRange, 0), (xRange, 0)], linestyle = '--', rgbcolor = (0, 0, 0))
    return c + L + R


def isInfinity(coordenadas):
    return coordenadas == [0, 1, 0]

def oposite(f, h, coords):
    if isInfinity(coords):
        return coords
    assert(len(coords) == 3)
    x = coords[0]
    y = coords[1]
    return [x, -y - h(x), coords[2]]

def isSpecial(f, h, point):
    coordenadas = point._coords
    return not isInfinity(coordenadas) and coordenadas == oposite(f, h, coordenadas)


def pointToString(point):
    if isInfinity(point._coords):
        return "infinito"
    return "(%d, %d)" % (point._coords[0], point._coords[1])

def pointsToString(pointsCurve):
    return map(pointToString, pointsCurve)

def specialPoints(pointsCurve, f, h):
    return filter(lambda point: isSpecial(f, h, point), pointsCurve)



def mumford(a, b, f, h):
    '''
    Returns a boolean indicating if [a, b] represents a divisor on the picard group 
    of the hyperellptic curve defined by f and h
    @param a: polynomial in one variable over a finite field
    @param b: polynomial in one variable over a finite field
    @param f: polynomial in one variable over a finite field
    @param h: polynomial in one variable over a finite field
    The curve is given by y^2 + h(x)y = f(x)
    '''
    genus = (f.degree() - 1) / 2
    assert a.degree() <= genus, "bad degrees: %d is greater than and %d" % (a.degree(), genus)
    semi_mumford(a, b, f, h)
    return True

def semi_mumford(a, b, f, h):
    '''
    Says if [a, b] represents a semi-reduced divisor on the picard group of the 
    hyperellptic curve defined by f and h
    '''
    assert a.divides(b ** 2 + b * h - f), "%r does not divide %r" % (a, b ** 2 + b * h - f)
    assert b.degree() < a.degree(), "bad degrees: a and b"
    return True

def semi_mumfordp(a, b, f, h):
    return a.divides(b ** 2 + b * h - f) and b.degree() < a.degree()

def irredutiveis(x, max):
    '''
    Prints all irreducibles of x^2 - a por a on the interval [0, max]
    '''
    for i in xrange(max):
        if (x ** 2 - i).is_irreducible():
            print "e irredutivel: x^2 - %d" % i

def addPointsTuples(p1, p2, f, h):
    return addPoints(p1[0], p1[1], p2[0], p2[1], f, h)

def addPoints(u1, v1, u2, v2, f, h):
    mumford(u1, v1, f, h)
    mumford(u2, v2, f, h)
    d1, e1, e2 = u1.xgcd(u2)
    d, c1, c2 = d1.xgcd(v1 + v2 + h)
    s1 = c1 * e1
    s2 = c1 * e2
    s3 = c2
    u = (u1 * u2) // (d ** 2)
    v = ((s1 * u1 * v2 + s2 * u2 * v1 + s3 * (v1 * v2 + f)) // d).mod(u)
    return (u, v)

def tryPoli(meup, generator, f, h):
    for poli in generator:
        try:
            if semi_mumford(meup, poli, f, h):
                print "poli is semi = %s" % poli
        except:
            pass

def searchSame(point, ringPolys, jacobian, f, h):
    for grau in xrange(0, 5):
        for poli in ringPolys.monics(grau):
            for graudois in xrange(0, grau):
                for polidois in ringPolys.monics(graudois):
                    if semi_mumfordp(poli, polidois, f, h) and jacobian((poli, polidois)) == point:
                        print (poli, polidois)


def listMumford(anelPolis, jacobiano, f, h):
    genus = (f.degree() - 1) / 2
    for grau in xrange(1, genus + 1):
        print "degree = %d" % grau
        for a in anelPolis.monics(grau):
            for b in polis(anelPolis, grau - 1):
                try:
                    mumford(a, b, f, h)
                    print "[%r, %r] is mumford of degree: %d" % (a, b, degreeM(jacobiano, a , b))
                except:
                    pass

def degreeM(jacobian, a , b):
    point = jacobian((a, b))
    for i in xrange(1, 15):
        if (point * i).is_zero():
            return i
    return - 1


def polis(polyRing, grau):
    return polyRing.polynomials(max_degree = grau)


def reduce(u, v, f, h):
    genus = (f.degree() - 1) / 2
    tu = (f - v * h - v ** 2) // u
    tv = (-h - v).mod(tu)
    while tu.degree() > genus:
        u = tu
        v = tv
        tu = (f - v * h - v ** 2) // u
        tv = (-h - v).mod(tu)
    return (tu.monic(), tv)
